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C D ' We propose a new method for simulating QCD at finite density, where interesting 

C " 3 ' phases such as the color superconductivity phase is conjectured to appear. The 

method is based on a general factorization property of distribution functions of 
■ observables, and it is therefore applicable to any system with a complex action. 

C D ' The so-called overlap problem is completely eliminated by the use of constrained 

simulations. We test this method in a Random Matrix Theory for finite density 
QCD, where we are able to reproduce the exact results for the quark number den- 
sity. The achieved system size is large enough to extract the thermodynamic limit. 
Q ' Our results provide a clear understanding of how the expected first order phase 

^ | transition is induced by the imaginary part of the action. We also discuss the non- 

, commutativity of the zero chemical potential limit (fi — > 0) and the thermodynamic 

limit, which is relevant to recent Monte Carlo studies at small u. 



, — i 

X 



1. Introduction 



Recently there are a lot of activities in QCD at finite density, where in- 
teresting phases such as a superconducting phase have been conjectured to 
appear 1 . At zero chemical potential Monte Carlo simulations of lattice 
QCD enable nonperturbative studies from first principles. It is clearly de- 
sirable to extend such an approach to finite density and explore the phase 
diagram of QCD in the T(temperature)-/i(chemical potential) plane. The 
main obstacle here is that the Euclidean action becomes complex once the 
chemical potential is switched on. 

Nevertheless QCD at finite density has been studied by various ap- 
proaches with exciting conjectures. First there are perturbative studies 
which are valid in the /i — > oo limit 2,3 . Refs. [4] and [5] use effective the- 
ories with instanton-induced four-fermi interactions. As for Monte Carlo 
studies two directions have been pursued so far. One is to modify the model 
so that the action becomes real. This includes changing the gauge group 
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from SU(3) to SU(2) 6 , and introducing a chemical potential with oppo- 
site signs for up and down quarks 7 . The other direction is to explore the 
large T and small fi regime of lattice QCD, where the imaginary part of 
the action is not very large 8 > 9 < 10 > 1:L . These studies already produced results 
relevant to heavy ion collision experiments, but more interesting physics 
will be uncovered if larger // regime becomes accessible by simulations. 

In Ref. [12] we have proposed a new method to simulate systems with a 
complex action, which utilizes a simple factorization property of distribu- 
tion functions of observables. Since the property holds quite generally, the 
approach can be applied to any system with a complex action. The most 
important virtue of the method is that it eliminates the so-called overlap 
problem, which occurs in the standard re-weighting method. Ultimately we 
hope that this method will enable us, among other things, to explore the 
phase diagram of QCD at finite baryon density. 

As a first step we test 1 the new approach in a Random Matrix Theory, 
which can be regarded as a schematic model for QCD at finite baryon 
density 14 . We also present preliminary results 15 , which reveal certain 
noncommutativity of the [i — > limit and the thermodynamic limit. 



2. Brute-force approach — reweighting method — 

Suppose we want to study the model defined by the partition function 

Z = j dUc- So+tr , (2.1) 

where So and T are real. Since the weight e~ So+ * r in (2.1) is not positive 
definite, we cannot regard it as a probability density. Hence it seems diffi- 
cult to apply the idea of standard Monte Carlo simulations, which reduces 
the problem of obtaining vacuum expectation values (VEVs) to that of tak- 
ing an average over an ensemble generated by the probability. One way to 
proceed is to apply the reweighting method and rewrite the VEV (O) as 

(Oe« r )„ 

(O) = ^ , (2.2) 

where the symbol ( • ) denotes a VEV with respect to the phase- quenched 
partition function 

Z = J dUc- So . (2.3) 

Since the system (2.3) has a positive definite weight, the VEV ( • ) can be 
evaluated by standard Monte Carlo simulations. However, the fluctuations 
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of the phase T in (2.2) grows linearly with the system size V. Due to huge 
cancellations, both the denominator and the numerator of the r.h.s. of (2.2) 
vanish as e ~ const - v as V increases, while the 'observables' e lT and Oe lT are 
of 0(1) for each configuration. As a result, the number of configurations 
required to obtain the VEVs with some fixed accuracy grows as e const y . 
This is the notorious 'complex-action problem'. Moreover when one sim- 
ulates the phase-quenched model (2.3), one cannot sample efficiently the 
configurations which are relevant to the calculation of the VEV (O). This 
is the so-called 'overlap problem'. 



3. New approach — factorization method — 

In the factorization method proposed in Ref. [12], the fundamental objects 
are the distribution functions (we assume the observable O to be real) 

p(x)^ f (S(x-O)) (3.1) 

p(°)(x) d ^ f (S(x-O)) (3.2) 

defined for the full model (2.1) and for the phase-quenched model (2.3), 
respectively. The important property of p(x) is that it factorizes as 

p(x) = ±pM(x)<p{x), (3.3) 

where the constant C is given by C A = (e* r )o- The 'weight factor' (p(x), 
which represents the effect of T, can be written as a VEV 

<p(x) = (e ir ) x (3.4) 
with respect to yet another partition function 

Z{x) = j dUc- Sa S(x - O) . (3.5) 

The 5-function represents a constraint on the system. In actual simulation 
we replace the (5-function by a sharply peaked potential. The distribution 
p^(x) for the phase-quenched model can also be obtained from the same 
simulation. Then the VEV (O) can be obtained by 

where the overlap problem is eliminated by forcing the simulation to sam- 
ple the important configurations by the constraint. The knowledge of the 
weight factor ip(x) is useful because it tells us precisely which values of O 
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are favored or disfavored by the effects of the oscillating phase. Once a 
rough estimate of p(x) is obtained, one may perform multi-canonical sim- 
ulations with an appropriate weight (instead of simulating (3.5) for many 
x) to sample relevant configurations more efficiently This has not yet been 
done, however. 



4. Random Matrix Theory for finite density QCD 

The Random Matrix Theory we study is defined by the partition function 



Z = J dWe~ Nt < w ^ dctD , 



(4.1) 



where W is a N x N complex matrix, and D is a 2N x 2N matrix given by 

m iW + 



D 



(4.2) 



The parameters m and ji correspond to the 'quark mass' and the 'chemical 
potential', respectively. The fermion determinant becomes complex for /i ^= 
0, so we write it as detD = e lF | detD|. The complex-action problem arises 
due to the phase T. In what follows we consider the massless case (m = 0) 
for simplicity and focus on the 'quark number density' defined by 
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(4.3) 



The model was first solved in the large- TV limit 14 , and turned out to be 
solvable later even for finite TV 16 . The partition function can be expressed 
as 

(_l)iv+i 



Z(p,N) = 7re K iV- (A,+1) AM 



1 + 



where k = —Nn 2 and 7(n,x) is the incomplete 7-function defined by 



j(n,x) = [ 
Jo 



From this one obtains the VEV of the quark number density as 
1 d 



= -M 



1 + 



n N e~ K 



{-l) N + 1 N+j(N + 1,k) 



(4.4) 

(4.5) 

(4.6) 
(4.7) 
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Taking the large- TV limit, one obtains 

lim (*) = {-? f r/i<Mc (4.8) 

where p, c is the solution to the equation l+^ 2 +ln(/i 2 ) = 0, and its numerical 
value is given by p c = 0.527- • • . We find that the quark number density 
(v) has a discontinuity at p, = p c . Thus the schematic model reproduces 
qualitatively the first order phase transition expected to occur in 'real' QCD 
at nonzero baryon density. 

The phase-quenched model defined by the partition function 

Z Q = J dWe- NtT( - w ' w ^\detD\ (4.9) 
can be solved in the large N limit 14 and one obtains 

IimMo = (? / !° rM< ! (4-10) 
a/->oo [ for fj, > 1 , 

which is a continuous function of the chemical potential p unlike (4.8). 
Thus the first order phase transition in the full model (4.1) occurs precisely 
due to the imaginary part T of the action. This model therefore provides a 
nice testing ground for simulation techniques for finite density QCD 14 . 



5. Testing the factorization method in the RMT 

Since v is complex for each configuration, we decompose it into the real and 
imaginary parts as v = z/r, + iv\ and calculate (z/r) and {v\) by the factor- 
ization method is purely imaginary). We introduce the distribution 
functions for z/r and v\ separately as 

Pi(x) =' (S(x - Vi)) i = R,l, (5.1) 
pr(*)= f (S(x-Vi)) i = R,I. (5.2) 
The factorization holds for both pn(x) and pi{x) as 

Pi{x) = ^ p\ 0) (x) <fii{x) i = R, I, (5.3) 

where the constant C is given by C = f (e ir )o- The weight factors tfi(x) 
can be written as a VEV 

<ft(z)= f (e ir )»,x i = R,I (5.4) 
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with respect to the constrained phase-quenched model 

Zi{x) = J dWe- NtI ^ w ^\detD\8{x-Ui) i = R, I . (5.5) 

Under the transformation W \— ► — W, the Gaussian action is invariant, 
whereas the fermion determinant det D as well as the quark number density 
v becomes complex conjugate. Due to this symmetry, we have 

!yS R (x)* = <fi R (x) , (5.6) 

= pi(-z) , (5.7) 
p { 1 \-x)=p[ 0) (x) . (5.8) 

Using these properties, we arrive at 

1 f°° 

{vr) = g J dx x p^ } (x) w R (x) , (5.9) 

(vi) = -£j o dxxp ( l 0) (x)w 1 (x) , (5.10) 

/OO 
dxp i ° ) (x)w R (x) , (5.11) 

where the weight factors Wi(x) are defined by 



w R (x) d = (cosr) R:X (5.12) 
wi(x) = {sinT)i. x = -wi(-x) . (5.13) 



Table 1. Results of the analysis of (v) described in the text. Sta- 
tistical errors computed by the jackknifc method arc shown. The last 
column represents the exact result for (u) at each /i and N . For /i = 0.2 
the exact result is (v) = —0.2 with an accuracy better than 1 part in 

10- 9 . 





N 


W> 


i(vi) 


<") 


(v) (exact) 


0.2 


8 


0.0056(6) 


-0.1970(5) 


-0.1915(7) 


-0.20000. . . 


0.2 


16 


0.0060(4) 


-0.1905(13) 


-0.1845(13) 


-0.20000. . . 


0.2 


24 


0.0076(9) 


-0.1972(14) 


-0.1896(17) 


-0.20000. . . 


0.2 


32 


0.0021(8) 


-0.1947(19) 


-0.1927(25) 


-0.20000. . . 


0.2 


48 


0.0086(37) 


-0.2086(54) 


-0.2000(88) 


-0.20000. . . 


1.0 


8 


0.8617(10) 


0.1981(13) 


1.0598(12) 


1.066501. . . 


1.0 


16 


0.8936(2) 


0.1353(6) 


1.0289(5) 


1.032240. . . 


1.0 


32 


0.9207(1) 


0.0945(2) 


1.0152(3) 


1.015871. . . 



In Table 1 we show our results for two values of fx, n = 0.2 and jj, = 1.0, 
which are on opposite sides of the first order phase transition point /i = 
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fi c = 0.527 • • • . They are in good agreement with the exact results, and the 
achieved values of N are large enough to extract the large N limit. 

In Fig. 1 we plot iur(j) for N — 8 at various /x. It is interesting that the 
wr(x) changes from positive to negative for \x < /i c , but it changes from 
negative to positive for it > /j c . (Similarly w\[x) is positive at x > for 
fi < He, but it is negative at x > for \i > /x c .) Thus the behavior of Wi{x) 
changes drastically as the chemical potential /x crosses its critical value /x c . 
These results provide a clear understanding of how the first order phase 
transition occurs due to the effects of T. Fig. 2 shows the results for (v) 
obtained by the factorization method for N = 8 at various /i, which nicely 
reproduce the gap developing at the critical point. 
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0.6139 
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0.8 - - - 
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1.0 - " 
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Figure 1. The weight factor u>b.(x) is plotted for N = 8 at various fi. 



6. Noncommutativity of /x — > and N — > oo 

In this Section we discuss the noncommutativity of the two limits, /j — ► 
and N — > oo. The absence of such noncommutativity is implicitly assumed 
in most of the recent approaches used in simulating finite density QCD at 
small /i. This includes the multi-parameter reweighting approach 8 , the 
Taylor expansion approach 9 , and the imaginary [i approach 10 - 11 . In all 
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Figure 2. The VEV (v) obtained by the factorization method is plotted against fi for 
N = 8 including the critical regime. Statistical errors computed by the jackknifc method 
are also shown. The dashed line represents the exact result (4.7) for (u) at N = 8. 

these works one restricts oneself to the small \x regime where the fluctuation 
of the phase is still under control. 

In fact the noncommutativity can be readily seen from the exact result 
(4.4) for the partition function (4.1). The phase of the determinant vanishes 
at /j, = for finite N, and one obtains a nonzero result for the partition 
function Z = 1 in the large N limit. On the other hand, the oscillation of 
the phase becomes pronounced at sufficiently large N even for small but 
finite fi, and as a result one obtains Z = in the large N limit as far as 
is kept finite. This implies in particular that the free energy 

/( M ) = -Jim ) ^lnZ( M ,iV) (6.1) 

has a discontinuity at \l — as 

lim i f(n) > /(0) = . (6.2) 

We expect that, in general, the free energy of a system with a complex 
action has a discontinuity at a point in the parameter space where the 
imaginary part of the action vanishes identically. 

With the factorization method we can take the two limits fj, — ► 0, 
N — > oo in different orders and compare the results. As we will see, we 
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do observe the noncommutativity in various ways. On the other hand, we 
know from the exact result (4.7) that the VEV (v) does not have the non- 
commutativity. In the factorization method (v) = (i/r) + i(v\) is calculated 
by the formulae (5.9), (5.10) and (5.11). In fact the functions w\(x) 
and (x) have the noncommutativity, but these effects cancel each other 
in the end results for (i/r) and {u\). In what follows we present preliminary 
results relevant to (vr), but similar results are obtained for (z/r) as well 15 . 




-1.5 -1 -0.5 0.5 1 1.5 



x 

Figure 3. The weight factor wr(x) is plotted for /i = 0.1 and 0.2 at N = 8, 16, 32. 



Let us first look at the weight factor Wr(x), which has the noncommu- 
tativity similar to the partition function. At fj, = one obtains wr(x) = 1 
for any N, whereas in the large ./V limit one obtains wr(x) = for any \x. 
In Fig. 3 we plot wr(x) for /i = 0.1 and p = 0.2 at TV = 8, 16, 32. It shows 
clearly that the behavior of m)r(i) depends much on the order of the two 
limits p — > 0, TV — ► co. 

Let us next turn to (x). In Fig. 4 we plot it for various N at fj, = 0.2. 
At small N the distribution is peaked near the origin and the dependence 
on TV is small. At sufficiently large N the peak moves to x ~ \i and starts 
to grow, which is consistent with the large result {isr)o = /i. Empirically 
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we find that the transition occurs at 

0.25-0.3 ,„„, 

N c = ^ ■ 6 -3 

M 

Thus the distribution p^\x) for the phase-quenched model also depends 
much on the order of the two limits \i — > 0, TV — > oo. This can be seen 
more clearly in Fig. 5, where we plot the VEV (^r)o against \x for various 
TV. In particular the derivative ^(i^r)o becomes if one takes the \i — > 
limit first, but it becomes 1 if one takes the TV — > oo limit first. 

The product p^\x)wr(x) gives the unnormalized distribution for vji 
in the full model, which we plot in Fig. 6. The distribution itself, even 
after appropriate normalization, has the noncommutativity, but the VEV 
(i/r) calculated by the formula (5.9) is always closed to zero (see Table 
1). The reason depends on the order of the two limits. If we take the 
TV — > oo limit first, the positive and negative regions of Pr(x) cancel each 
other in the calculation of the first moment. If we consider small \i first, 
the distribution is peaked around the origin, which makes the first moment 
close to zero. Thus the noncommutativity cancels in the end result for the 
VEV (i/r). 
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Figure 4. The function pj^'(x) is plotted for fi = 0.2 at various N. 
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7. Concluding remarks 

The factorization method has been applied also to other systems with com- 
plex actions. In the original paper 12 , it was used to study the dynamical 
generation of space time in superstring theory based on its matrix model 
formulation 17 . In this case the weight factors turned out to be positive 
definite, which enabled us to use their scaling property to make an ex- 
trapolation to larger system size. The method 18 proposed for simulating 
9- vacuum like systems can be regarded as a special case of the factorization 
method. Promising results are obtained in the 2d CP 3 model etc. 
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